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ABSTRACT 


The nonlinear temporal evolution of disturbances in compressible flow be- 
tween infinitely long, concentric cylinders is investigated through direct nu- 
merical simulations of the full, three-dimensional Navier-Stokes and energy 
equations. Counter-rotating cylinders separated by wide gaps are considered 
with supersonic velocities of the inner cylinder. Initially, the primary distur- 
bance grows exponentially in accordance with linear stability theory. As the 
disturbances evolve, higher harmonics and subharmonics are generated in a 
cascading order eventually reaching a saturation state. Subsequent highly 
nonlinear stages of the evolution axe governed by the interaction of the dis- 
turbance modes, particularly the axial subharmonics. Nonlinear evolution of 
the disturbance field is characterized by the formation of high-shear layers 
extending from the inner cylinder towards the center of the gap in the form of 
jets similar to the ejection events in transitional and turbulent wall-bounded 
shear flows. 
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1 INTRODUCTION 


The interaction of viscous and centrifugal mechanisms in the stability of high- 
speed shear flows can be studied in a relatively simple model offered by the 
compressible rotating Couette flow. In the present work, direct numerical 
simulations are performed to investigate the nonlinear temporal evolution of 
forced disturbances in wide gaps between counter-rotating, concentric cylin- 
ders of infinite length. 

The balance of viscous and centrifugal forces determines the first transi- 
tion in the incompressible Couette flow which can be predicted by the linear 
stability theory . 1 Further interaction of these two mechanisms generates rich 
complex flow patterns 2 which render the problem attractive but the ensu- 
ing nonlinearities are hardly amenable to a thorough theoretical treatment. 
According to Stuart 3 nonlinear effects can be the generation of harmonics 
(and subharmonics) of the fundamental mode, distortion of the mean motion, 
and the distortion and evolution of the fundamental mode as well as of the 
other modes. Finite-amplitude incompressible Taylor- Couette flow has been 
treated in the vicinity of the neutral-stability region by Davey 4 by weakly- 
nonlinear methods based on the amplitude-expansion procedures suggested 
by Stuart . 5 One critical assumption used by Davey is that all harmonics of 
axial periodicity add in phase. Snyder and Lambert 6 demonstrated in their 
experiments that all harmonics are indeed in phase with the fundamental 
resulting in the emergence of jets in the flow field as a direct consequence of 
nonlinear interactions. Nonlinear jet-modes were also identified by Lorenzen, 
Pfister & Mullin 7 in the experiments of incompressible Taylor- Couette flow 
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with a stationary outer cylinder. 

Weakly nonlinear theories consider a subset of nonlinear effects, include a 
limited number of modes and are often restricted to the proximity of the re- 
gion where the disturbances are only marginally unstable. Among such stud- 
ies, Jones 8,9 used a mixed finite-difference/spectral technique to study the 
nonlinear incompressible Taylor-Couette flow effectively for Reynolds num- 
bers (or Taylor numbers) well beyond the critical values. He found the jet- 
mode to be the most unstable nonlinear mode as opposed to the classical 
wavy mode for some parameter ranges. Jones also found two distinct modes, 
i.e. the harmonic and subharmonic jet-modes. The subharmonic mode is 
slightly more unstable; therefore the harmonic mode has not been detected 
in experiments. 

The study of compressible flow instabilities is more complicated than the 
incompressible ones because disturbances in density, temperature, viscosity 
and conductivity have to be considered in addition to pressure and velocity 
fluctuations . 10 Hatay et al. n studied the linear stability of compressible 
Taylor-Couette flow at finite Mach numbers for three-dimensional distur- 
bances. Accordingly, compressibility can either stabilize or destabilize the 
growth of small-amplitude disturbances depending on the gap width and the 
rotational speeds of the cylinders. By solving the fully nonlinear, compress- 
ible equations, the present study aims at extending the previous results from 
the linear theory including all possible nonlinear effects in a parameter range 
where the disturbances are predicted to be linearly most unstable. 
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2 GOVERNING EQUATIONS 


We consider viscous compressible flows governed by the Navier-Stokes and 
energy equations and the thermodynamic equation of state. These equations 
are written for a heat conducting perfect gas in a three-dimensional cylindri- 
cal coordinate system, ( r,0,z ) denoting radial, azimuthal and axial directions 
respectively. The flow field is bounded by two infinitely-long, concentric 
cylinders; we define R as the radius and Cl as the angular velocity of the 
cylinders (Figure 1). Subscript ()i refers to the inner cylinder and subscript 
() 2 to the outer cylinder; superscript ()* is used for dimensional quantities. 
Scalar flow variables are density p, temperature T, and pressure p and the 
velocity vector is denoted by u = (tt r , ug, u 2 ) . Length, time, velocity, density, 
pressure and temperature are scaled with Rf, ftp 1 , ft* JEJ, p*, p'(R^Cll) 2 and 
Tj* respectively. Thermophysical properties, viscosity and heat conductivity 
(p* and k*), are also scaled with their respective values on the surface of the 
inner cylinder. 

The governing equations can be expressed in terms of the dimensionless 
variables in vector form as follows: 

^ + V-(pa) = 0, (1) 

^ + V.(„uu) = -V P + iv.., (2) 

£|I + V • (J5ru) + ■ (,VT) 

= Te V ' ^ ~ V ‘ ^ 

*yM 2 p = pT, 
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( 3 ) 

( 4 ) 


where 


r = /i(V u + Vu T ) + A(V • u)I. 


( 5 ) 


is the viscous stress tensor and A = —2/3/i is the second coefficient of viscos- 
ity; () T is the transpose operator, and I is the identity tensor. Total energy 
is the sum of kinetic and internal energies: 
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Reynolds number, Mach number and Prandtl number above are defined as: 


Re = 


i*i/ pi ' 


M = 


•J'iR-T;' 


Pr = 


PjCp 

ki 


( 7 ) 


where R* = c* — c* is the gas constant. Variation of Prandtl number Pr, and 
specific heats c*, c* with temperature is neglected. In this study Pr and spe- 
cific heat ratio 7 =c*/c* are taken to be 0.7 and 1.4 respectively. Sutherland’s 
law is used to express the dependency of viscosity /i on temperature, 


P 


= k = T* 


1 + 5 
r + 5’ 


(8) 


where 5 = 110AK/T^ and T* = 288.15tf. 

No-slip velocities and prescribed temperatures are imposed along the ra- 
dial direction on the cylinder surfaces. The boundary conditions for the 
pressure are obtained from the radial momentum equation. Periodicity is 
assumed in the azimuthal, 6, and the axial, 2 , directions. The assumption 
of axial-periodicity is equivalent to following the evolution of disturbances 
in time with respect to a reference frame moving with the phase (or group) 
velocity of the disturbances. Consequently, the computational box-length 
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corresponds to the largest wavelength represented in the solution. The choice 
of the box-length is not a trivial task as the most critical waves should be 
included in the simulation. For this purpose, results from an earlier linear 
stability analysis of this problem 11 were used to determine the wavelength 
of the most critical disturbances. 

The base flow is a steady solution of the governing equations. For this 
purpose, axisymmetric two-dimensional flow in the (r,0) plane is assumed 
so that all variables depend only on the radial coordinate, r. Velocity com- 
ponents in the axial direction, z, and the radial direction, r, are taken to 
be zero for a purely azimuthal base flow. The resulting system of ordinary 
differential equations are then solved as a boundary value problem. 11 

The general form of the initial distributions of the flow variables used in 
the simulations is: 


$(r,0,z,t = 0) = $ 0 + + A 2 $2 + (9) 

where $ = (p, pu T , pug, pu 2 , Et) t is the vector of primary variables, and the 
base flow is denoted by $o- The primary disturbance, is obtained from the 
least-damped linear, complex eigensolution, with a maximum amplitude 
of Ai: 

= Real [$i(r) exp[i(az + /?0)] j , (10) 

where a and /? are the wavenumbers in the z and 6 directions, respectively, 
and i = \J— 1. The axial subharmonic of the primary disturbance is similarly 
obtained from the linear stability analysis as $2 with the wavenumber pair 
(a/2,/?) and is introduced with a maximum amplitude of A 2 into the initial 
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field: 


$ 2 = Real $ 2 (r) exp (i(^ z + Z^))] , (11) 

Note that zero phase is used for $1 and $ 2 in Equations 10 and 11. A 
random field, at a maximum amplitude of Ar is also added to the ini- 
tial disturbance fields to trigger the nonlinear evolution of the instabilities. 
The random field is generated in the physical space by using random num- 
ber generators provided by the Cray Scientific and Mathematical Functions 
Library. 

The computational domain length in the axial direction, L z , is specified 
as 2n /a in the harmonic simulations and as 4ir/a when the evolution of 
axial-subharmonics is investigated. 

3 NUMERICAL METHOD 


The solution procedure used in the present work is the Two-Four method 
proposed by Gottlieb & Turkel 12 as a fourth-order variant of the fully-explicit 
MacCormack method. The Two-Four scheme is highly phase- accurate and is 
very suitable for wave propagation and wave interaction problems. 13 In the 
present work, the Two- Four method is adapted to the cylindrical geometry 
with a source-like treatment of the centrifugal terms. A stretched mesh in 
the r direction is used to cluster more nodes towards the solid walls. The 
allowable time step, A t, for the current fully explicit method depends on the 
Courant Number ( CFL ) which is defined as: 


CFL = At 


max 

grid 


.^li I JM . M _L a I 1 I 1 I 1 

Ar rAO A z \J Ar 2 (r A0) 2 Az 2 
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( 12 ) 


+ 


0 vis ( 1 1 1 \ 

Rt + (rA0) 2 + A?JJ ’ 


where the term vis is given as, 


max 


vis = 


(/i,A + 2p, 


and a is the local speed of sound, 


a = 



(13) 


(14) 


The stability of the numerical scheme is ensured for CFL < 1 and the scheme 
is optimized for CFL « 2/3. 

The solution procedure was implemented on the CRAY-YMP at NASA 
Langley Research Center taking full advantage of the vector architecture of 
this processor. It is desired that in the absence of disturbances, the base flow 
should satisfy the discretized full governing equations to the order of machine 
accuracy available by the computing platform. In the present problem, this 
is achieved by introducing extra source terms (which represent the base flow) 
into the discrete form of the 0-momentum and total energy, Ej, equations. 
These source terms are calculated simply as the difference between the left 
and right-hand sides of the discrete equations (typically of order 10 -12 ) and 
axe kept constant thereafter throughout the computation. This technique 
was applied by Erlebacher & Hussaini 14 in their simulations of supersonic 
boundary-layer flow. For the present flow field, this operation is equivalent 
of obtaining fully nonlinear, complete perturbation equations, hence the be- 
havior of the system with or without the source terms does not change in 
any way. 
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4 


RESULTS 


4.1 Initial Disturbances 

The Navier-Stokes solver described in Section 4 was tested extensively by 
comparisons with the linear theory; these results are reported elsewhere. 11 
Accordingly, the amplitude and phase agreement between the linear theory 
and the direct numerical simulation is better than 1.0%. 

In the present work, we investigate the nonlinear temporal evolution of 
compressible flow in a wide gap (/?i/f? 2 = 0.5) between infinitely-long, con- 
centric cylinders with equal and opposite rotation speeds (Q 2 /fIi = — 1). 
Direct numerical simulations were performed at a supersonic Mach number 
of M = 2.0 and at a Reynolds number of Re = 2, 000. According to the linear 
stability analysis 11 the primary instability is a three-dimensional, traveling 
wave when the cylinders are in the counter-rotating configuration at finite 
Mach numbers. At Re = 2,000 a wide band of linearly unstable distur- 
bances exist with different axial and azimuthal wavenumbers (a, j3) and the 
combination of a = 19.0 and — 1 gives the highest growth rate; due to 
the circular geometry, f3 can assume only integer values. For these parameter 
values, a high-resolution linear analysis with 121 nodes in the radial direction 
generates the temporal eigenvalue (complex frequency) u = 0.5168 + *0. 9758, 
whose imaginary part is the growth rate. The corresponding eigenfunctions 
(Figure 2) were used as the primary disturbance field ($i in Eq. 11) for the 
Navier-Stokes simulations and reveal a one-cell structure with maxima closer 
to the inner cylinder and quiescent behavior in the outer half of the domain. 
The eigenfunctions for the axial subharmonic ($ 2 in Eq. 11) were obtained 
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in a similar manner and are presented in Figure 3. 

In the subsequent sections, we report the results from three separate di- 
rect numerical simulations. The first simulation discloses the evolution of the 
harmonics in the absence of the subharmonic generation, whereas the second 
allows for the evolution of the subharmonics. The last simulation is initial- 
ized with comparable finite amplitudes of the primary and axial-subharmonic 
disturbances with zero phase difference, revealing the effects of forced sub- 
harmonics. We analyze the nonlinear development of the disturbance field 
for each case. 

4.2 Harmonic Development 

This simulation was initialized with the primary disturbance imposed on the 
base flow with a maximum amplitude of A\ = 10~ 3 together with a random 
disturbance field of amplitude Ar = 10 -6 . The computational domain was 
chosen to accommodate one full axial wavelength of the fundamental distur- 
bance ( L z = 2tt / a). The temporal evolution of the disturbances was then 
studied with a computational grid resolution of 65 x 66 x 34 in the radial, 
azimuthal, and axial directions, respectively. During the course of the sim- 
ulation, the solution fields were interpolated onto a high-resolution grid of 
65 x 66 x 66 and integrated in time at this resolution. The complete agreement 
between these two simulations indicate that solutions are grid independent. 

In order to analyze the data, we first performed one- dimensional Fourier 
expansions in one periodic direction (either 6 or z) and followed the time 
evolution of the maximum Fourier amplitudes. For example, an expansion 
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of the field variable <f> in a one-dimensional Fourier series along the axial 
direction, z, can be written as: 

<t>( r ,z,0,t) = exp (»*!<«), (15) 

*1 

where |&i| < \N Z and N z is the number of nodes in the z direction. The 
maximum amplitude is found in the (r, 9) plane: 

= cm, (i6) 

where ()* denotes the complex conjugate. Similarly, a one-dimensional Fourier 
expansion in the azimuthal direction , 0, yields: 

<t>(r,z,9,t) = 5Z k r , z, k 2 , t) exp (ik 2 /36), (17) 

h 

where \k 2 \ < \Ng and N g is the number of nodes in the 9 direction. This 
time, the maximum amplitude is found in the (r, z) plane: 

(#*). (is) 

As predicted by the linear theory, the initial growth of the primary dis- 
turbance is exponential. In Figure 4, the time evolution of A(ki = l,t) and 
A(k 2 = 1, t) is compared with the linear theory. The quantitative comparison 
of the growth rates in Table 1 reveals good agreement in the initial stage. 

The time history of A(k\,t) and A(k 2 ,t) is presented in Figure 5 for the 
total energy, Ej. In each of these plots, the topmost curve represents the 
first wavenumber {k\ = 1 or k 2 = 1 in Eq. (16) ) and the lower curves are 
the harmonics in ascending order. The amplitudes of the harmonics first 
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decay and then are energized and follow the rapid growth of the primary 
mode. The excitation of the harmonics takes place in a perfect hierarchy in 
which the smallest wavenumber becomes energized first, reaches saturation 
at a finite amplitude and then transfers its excess energy to the next smaller 
wavenumber. This cascading process continues until all the modes reach a 
nonlinearly quasi-saturated plateau stage along which the disturbances are 
sustained at finite amplitudes. The whole range of harmonics, not a few 
wavenumbers, in both the axial and azimuthal directions is supported in this 
stage. Power spectra of A, the maximum amplitude of the total energy Et 
as defined in Eqs. (16) and (17), at t = 25.51 (Figure 6) are smooth and 
continous with a fast decay towards the large wavenumbers indicative of the 
sufficient grid resolution used in the numerical solution at this stage. The 
smoothness of the spectra is related to the continous cascade of energy from 
the smaller wavenumbers. 

Figure 7 illustrates the generation of the harmonics in the two-dimensional 
wavenumber space. Here, Et is expanded in a Fourier series in 0 and z di- 
rections as: 

E T (r,z,6,t) = ^T{r, h, k 2 ,t) exp (ikiaz + ik 2 00), (19) 

ki k 2 

where |fci| < \N Z , |fc 2 | < \Ng, N z and Ng are the number of nodes in z and 
9 directions, respectively. The amplitude of these modes are defined as: 

A(ki,k 2 ,t) = — — f EtEj dr. (20) 

fi 2 — Hi 

In the present analysis, the quantity given by Eq. (20) refers to the energy 
contained in the spectral mode (ki,k 2 ). This definition of modal energy is not 
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unique and differs from the incompressible definition. However, the present 
approach has the advantage of containing both the kinetic and internal en- 
ergies of the disturbance fields. Moreover, the analysis based on the present 
definition gives a coherent and consistent picture of the events in the physical 
space as well as in the spectral space. 

As depicted in Figure 5, the harmonics are energized sequentially as time 
progresses. Also, new peaks appear in the two-dimensional Fourier spectrum 
shown in Figure 7 in a systematic fashion along the axial-harmonics. The 
direction of this hierarchial generation of harmonics is initially: 

(na,n/3) — ► [(n + l)a, (n + l)/3] (21) 

where n = 1, 2, ...\N Z . Once the harmonics are energized in this hierar- 
chy, energy is transferred to the other wavenumbers in the two-dimensional 
wavenumber space. At later times, a uniform distribution of the amplitudes 
is approached with gradual decay towards the smaller wavenumbers. It is 
interesting to note that the amplitude-decay is much slower in the azimuthal 
direction than in the axial direction. Therefore, the azimuthal mesh resolu- 
tion requirements became prohibitive beyond t = 40.0. 

The main energy source in this temporal evolution is the external torque 
applied at the rotating cylinders. Figure 8 presents the mean velocities (av- 
eraged in the z-9 plane) at selected times. The increase in the gradients close 
to the i nn er cylinder is due to more efficient energy and momentum transfer. 
Relatively flat portions of the ug profiles for about r « 1.1 indicate enhanced 
momentum- mixing in that region. The development of non-zero u r and u z 
is related to the nonlinear distortions of the initially sinusoidal disturbance 
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fields. Consequently, in the inner half of the gap, the mean flow is no longer 
purely azimuthal but has a weakly swirling component. 

The events in the spectral and physical spaces can be interconnected 
through the inspection of instantaneous disturbance fields. Figure 9 presents 
contour plots of disturbance total energy, Et, in several (r-z) planes. During 
the early stage, t = 0.6313, the initially imposed randomness quickly dies 
and the fundamental mode with its characteristic sinusoidal variation dom- 
inates the flow field at this cross-section. At later times, as the harmonics 
are energized and the plateau stage is reached, the formation of a jet-like 
structure is observed close to the inner wall. The harmonics are generated in 
phase with each other; consequently, their superposition leads to this high- 
shear-layer formation. The jet-like structure expands towards the center of 
the gap in the shape of a mushroom but does not quite extend into the 
outer half of the gap so that this region remains stable. Along the plateau 
stage ( t > 20), these jets or mushroom-like ejection structures promote even 
higher-shear layers with stronger gradients close to the inner wall. Three- 
dimensional isosurfaces of disturbance total energy, Et, at the levels of 50% 
of the maximum of the disturbance field, are presented in Figure 10 to illus- 
trate this sequence of events in greater detail. In the beginning of the plateau 
(t = 12.85), the isosurfaces are perfectly spiral, indicating the dominance of 
the fundamental mode. The development of high shear layers is discernible in 
the high energy isosurface as it narrows close the inner cylinder and spreads 
towards the center of the channel. Along the plateau (t = 31.72), the jets 
get thinner and the structures become susceptible to the higher instability 
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mechanisms as evident from the waviness superimposed on the spirals. This 
critical wavy development in the azimuthal direction is also evident from 
the increasing amplitudes and the crowdening if the azimuthal spectrum for 
t > 30 in Figure 5. Eventually (t=35.48), the spiral symmetry is lost as the 
high energy structure locally detaches from the wall generating small-scale 
structures during the breakdown process. The investigation of the physics 
of flow transition beyond this stage would require an extremely high spatial 
and temporal resolution not affordable by the computer resources available 
for this study. 

4.3 Subharmonic Development 

In this simulation, the computational domain is extended to accommodate 
two axial- wavelengths of the fundamental disturbance ( L = 2\ z — 4ir/a) in 
order to allow the generation and the evolution of the axial-subharmonics of 
the fundamental mode. For this case, the number of nodes is increased to 
81 x 130 x 102 in the radial, azimuthal, and axial directions, respectively. 
The other parameters and the initial conditions are kept the same as in the 
previous harmonic simulation. The first subharmonic mode is excluded from 
the initial conditions (A 2 = 0) and can be generated only through nonlinear 
mechanisms. 

Figure 11 depicts a similar evolution of the amplitudes of the one-dimensional 
Fourier components as in the harmonic case: initial rapid growth of the fun- 
damental mode saturates to reach a plateau stage. For all the simulations 
performed in the present study, the nonlinear saturation is achieved over ap- 
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proximately one rotation of the cylinders (period = 27r/fl 1 = 6.28). There- 
fore, the mechanism for the generation and the immediate evolution of both 
harmonic or subharmonic disturbances must be the inviscid, rotational insta- 
bilities. As the disturbances evolve and attain finite amplitudes, a nonlinear 
saturation state is reached. 

The evolution of the subharmonics is distinct from that of the haxmonics. 
Starting from the small-amplitude random initial disturbance field, subhar- 
monic modes neither reach the saturation level nor do they catch up with 
the harmonics until the very late stages of the transition process. 

Figure 12 compares the spectral time evolution of the axial-harmonic 
modes from a preliminary run performed up to t = 21 at a resolution of 
65 x 66 x 102 with the 81 x 66 x 102 high resolution case. This comparison 
reveals a similar behavior providing further justification to the adequacy of 
the computational grid used throughout the present study. 

Figure 13 shows the radial distributions of the Fourier amplitudes for 
the mean flow, the fundamental mode, and the first subharmonic mode at 
different times. Maximum amplitudes occur close to the inner cylinder, and 
strong development of a secondary extremum for r ~ 1.2 is observed at later 
times. The amplitudes at the inner cylinder surface are not zero, because as 
the kinetic energy of the disturbances diminishes to zero, the internal energy 
content actually increases very rapidly in the near-wall region. 

In Figure 14, the subharmonic components attain a well-developed spec- 
trum showing the decay of energy towards higher wavenumbers. The spectra 
of subharmonics and harmonics are similar but separate, suggesting a dif- 
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ferent route for energy transfer for the subharmonics. A nonlinear energy 
transfer analysis can usually be instrumental in clarifying interactions be- 
tween different modes 15 and can be performed by expanding the terms of 
the energy equation (Eq. 3) in Fourier series (Eq. 19) for the temporal 
evolution of the present flow. Neglecting the viscous terms because their 
contribution to mode interactions is secondary to the convective terms, the 
following equation is obtained from Eq. (19) for the nonlinear evolution of 
the Fourier modes: 

r ’£’ — - = - E E fr(r, k 2 , k[,k'v *) 

Ot . / u 

*1 K 2 

where T is the energy transfer rate to the mode (ki,k 2 ) from {k[, k' 2 ) and 
(ki — k[, k 2 — k 2 ) modes: 

T(r,ki,k 2 ,k[,k! 2 ,t) = -V • [#r(r, k[, k' 2 , t)u(r, h - k[,k 2 - k' 2 ,t) 

+p(r,k[,k' 2 ,t)u(r,ki - k[,k 2 -k 2 ,t)]. (23) 

The interaction between the modes can now be quantified as: 

T(k u k 2 , KX,t)= T(T)ir. (24) 

Note that this is a positive-definite quantity and does not show the direc- 
tion of transfer. The interactions of the harmonic, ( na,n 0), and the axial- 
subharmonic, ([n — l/2]a, n(3), modes with the mean flow, mode (0,0), is 
presented in Figure 15. The evident similarity between energy transfer rates 
displayed in these figures and the mode amplitudes shown in Figure 14 reveals 
that there is a direct energy transfer to each mode from the mean flow. When 



dE T { 


18 


a harmonic mode reaches saturation so does its interaction with the mean 
flow. Along the saturation state, the mean flow continues to supply energy to 
the modes and the modes which reach saturation transfer their excess energy 
to the higher modes down the nonlinear energy cascade. In the meantime, 
the subharmonic modes never reach saturation and their interaction with the 
mean flow continously increases during the simulation. 

Interactions of the fundamental mode, (a, /?), with the other modes (Fig- 
ure 16) indicate a similar energy transfer mechanism for the fundamental 
compared with the mean mode. Accordingly, the fundamental mode consti- 
tutes the next level in the harmonic cascade after the mean flow in trans- 
ferring energy to the harmonic and subharmonic modes and this scenario 
repeats itself along the harmonic cascade. Considering all modal energy 
transfer combinations, when a mode reaches saturation, net energy transfer 
to that mode diminishes to zero as the energy cascade continues to operate 
in equilibrium. 

Figure 17 displays the time evolution of the disturbance total energy, Ej, 
in several (r, z) sections for 9 = 0 and 7 r. Following the initial exponential- 
growth stage, jet-like ejection structures develop as in the harmonic case. 
These structures occur with the periodicity of the fundamental disturbance 
and are indeed harmonic in nature. At later times, subharmonic structures 
develop and compete with the harmonics. The emergence of subharmonic 
structures in physical space is the result of merging adjacent ejection struc- 
tures yielding a periodicity with twice the wavelength of the fundamental 
disturbance in the axial direction. These subharmonic jet-like structures 
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are stronger than their harmonic predecessors and reveal a powerful ejection 
extending from the wall region spreading towards the outer cylinder. The 
subharmonic structures quickly become susceptible to the higher instabili- 
ties leading to flow breakdown and are also accompanied by changes in the 
temporal response. In Figure 18, time-history data collected at a fixed point 
in space are presented. Quasi-periodic behavior along the plateau stage is 
discernible by two sequential peaks followed by two troughs with a period 
close to that of the fundamental mode (2t / R eal(u) = 12.16). The distur- 
bances have gained appreciable amplitude of about 0.05-0.10 at this stage. 
This temporal response is nonlinearly distorted forming spikes concurrent 
with the formation of jet-like structures in space. The loss of spatial sym- 
metry takes place simultaneously with the loss of quasi-periodicity in time, 
indicating the end of the plateau stage. 

The near-wall events investigated throughout the present study can also 
be traced from the the development of equishear lines at the wall. Figure 19 
illustrates the development of the disturbance part of fidug/dr at the inner 
cylinder wall. Initially, the equishear lines are straight with an inclination of 
a//? in the ( z-9 ) plane tracing the spiral structure imposed by the fundamen- 
tal disturbance. Along the plateau stage, the formation of high-shear layers 
caused by the strong jet-like structures is indicated from the intensity of the 
equishear lines. The flow becomes susceptible to higher instability mecha- 
nisms evident from the waviness of these contours at later times and the 
breakdown follows the emergence of structures with higher wavenumbers. 
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4.4 Forced Subharmonic Development 

The simulation described in Section 4.3 indicated that when subharmonics 
are allowed to develop in the computational box by prescribing the appro- 
priate box-length, they may dominate over the harmonic structures. In the 
previous simulations, although linearly the most unstable mode was chosen 
as the primary disturbance, the axial- subharmonics clearly have a crucial 
role during the nonlinear stage. To further clarify this issue, in the third 
simulation we introduced the subharmonic disturbances in the initial condi- 
tions with A 2 = 10~ 3 together with the harmonic disturbances ( Ai = 10~ 3 ); 
a random field ( Ar = 10~ 6 ) was also imposed on the base flow. 

Starting with equal amplitudes for the first subharmonic mode and the 
fundamental, during the exponential-growth stage (Figure 20) the axial- 
subharmonic modes evolve in a manner similar to the evolution of the har- 
monics. Both the subharmonic and the harmonic modes reach finite-amplitude 
levels for nonlinear saturation after approximately one rotation of the cylin- 
ders. However, the plateau stage is much less well-defined and less stable than 
the previous simulations and amplitudes of the azimuthal modes fluctuate 
along the plateau. In Figure 21 the amplitude spectra of the one-dimensional 
Fourier modes indicate slower decay towards the small wavenumbers com- 
pared with the previous simulations and the spectrum of the azimuthal modes 
is no longer smooth but contains small-scale oscillations. Consequently, the 
forced subharmonics appear to have a detrimental effect on the sustainability 
of the plateau stage. 

As the amplitudes of the harmonics and subharmonics become compara- 
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ble, nonlinear developments lead directly to the formation of jet-like struc- 
tures of the subharmonic type; no initial harmonic stage is observed during 
this simulation. Moreover, these structures are no longer confined to the 
inner half of the gap but extend further into the outer half which remained 
stable during the nonlinear evolution in the previous numerical simulations 
(Sections 4.2 and 4.3). Although the jets do not quite impinge upon the 
outer cylinder, they are strong enough to instigate significant activity in the 
outer half of the gap (Figure 22). 

As in the previous case, the formation of high-shear layers is discernible 
from the intensification of equishear lines in the contours of the disturbance 
part of [idue/dr at the inner cylinder wall (Figure 23). Characteristic of this 
simulation is the increasing subharmonic character of the flow field during 
the plateau stage. The resulting structures quickly become subject to higher 
instability mechanisms and lose their spiral symmetry. 

In all of the cases presented in this work, sharp boundary layers form 
on both the outer and inner cylinder walls during the nonlinear stage; how- 
ever, the ejection type events originate only close to the inner cylinder giving 
further indication that they are produced mainly by rotational instability 
mechanisms. As the nonlinear evolution progresses, rotational instabilities 
interact and compete with viscous instabilities. Instabilities occuring in the 
outer region in the forced-subharmonic case is likely to be a product of such 
interactions. In all stages of the evolution, compressibility maintains a reces- 
sive role and does not cause any peculiar changes in the nature of instabilities 
which are nurtured by viscous and centrifugal mechanisms. 
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5 CONCLUSIONS 


Three-dimensional direct numerical simulations of the nonlinear evolution 
of disturbances in compressible rotating Couette flow were performed at a 
moderate Mach number (M = 2.0) and at a Reynolds number of Re = 
2000.0. Initial conditions consisted of a random field imposed on the small- 
amplitude fundamental and the first axial-subharmonic modes which were 
obtained as three-dimensional traveling wave solutions from a high-resolution 
linear stability analysis. 

The initial growth of the fundamental mode obtained from the present 
DNS is in good agreement with the linear stability predictions. The har- 
monics and subharmonics first decay in favor of the exponentially-growing 
fundamental mode and then become excited and start to grow in the or- 
der of increasing wavenumbers. The growth of the modes saturates when 
finite amplitudes are reached; this nonlinear saturation state is referred to 
as the plateau stage. When excluded from the initial conditions, the axial- 
subharmonics do not attain finite amplitudes to reach the plateau stage. For 
initial conditions that contain equal amplitude forced subharmonic and the 
fundamental disturbances, the subharmonics reach a nonlinearly saturated 
state; however, this plateau stage cannot be sustained. Our results suggest 
that the axial subharmonic development is the path of transition dominating 
the nonlinear evolution of this flow. 

Along the plateau stage jet-like structures form carrying the fluid from the 
inner cylinder to the center of the gap. In the absence of forced subharmonics, 
these ejection events are harmonic structures and tend to merge to attain a 
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subharmonic character as they evolve. Otherwise, for forced subharmonics, 
subharmonic jets dominate both the inital and later stages of the transition 
process. The jets grow stronger and narrower and become susceptible to 
higher instability mechanisms. Eventually, spiral symmetries are broken by 
the detachment of the jets from the wall in a non-uniform manner generating 
small-scale structures during the process. The end of the plateau stage is also 
marked by the loss of quasi-periodicity in the temporal response. 

The events during the transition process of this supersonic shear flow are 
mainly generated by the interaction of the centrifugal and viscous mecha- 
nisms which are prevalent in incompressible shear flows as well. Only in- 
directly does compressibility affect the essentially incompressible transition 
phenomena in this shear-driven flow. 

References 

1 G. I. Taylor, “Stability of a viscous fluid contained between two rotating 
cylinders”, Phil. Trans. A 223, 289 (1923). 

2 D. Coles, “Transition in circular Couette flow”, J. Fluid Mech. 21, 
385 (1965). 

3 S. T. Stuart, “On the nonlinear mechanics of wave disturbances in stable 
and unstable parallel flows, part 1. the basic behavior in plane poiseuille 
flow”, J. Fluid Mech. 9, 353 (1960). 

4 A. Davey, “The growth of Taylor vortices in flow between rotating cylin- 
ders”, J. Fluid Mech. 14, 336 (1962). 


24 


5 S. T. Stuart, “On the non-linear mechanics of hydrodynamic stability”, 
J. Fluid Mech. 4, 1 (1958). 

6 H. A. Snyder and R. B. Lambert, “Harmonic generation in Taylor vortices 
between rotating cylinders”, J. Fluid Mech. 26, 545 (1966). 

7 A. Lorenzen, G. Pfister, and T. Mullin, “End effects on the transition 
to time-dependent motion in the Taylor experiment”, Phys. Fluids 26, 
10 (1982). 

8 C. A. Jones, “On flow between counter-rotating cylinders”, J. Fluid Mech. 
112, 443 (1982). 

9 C. A. Jones, “The transition to wavy Taylor vortices”, J. Fluid Mech. 
157, 135 (1985). 

10 L. M. Mack. “Boundary-layer linear stability theory”, Special Course on 
Stability and Transition of Laminar Flow, edited by R. Michel (AGARD, 
Neuilly sur Seine, France, 1984), AGARD Report No. 709, p. 3-1. 

11 F. F. Hatay, S. Biringen, G. Erlebacher, and W. E. Zorumski, “Stability 
of high speed compressible, rotating Couette flow”, Phys. Fluids A, 5, 
393 (1993). 

12 D. Gottlieb and E. Turkel, “Dissipative two-four methods for time- 
dependent problems”, Math. Comput. 30, 703 (1976). 

13 S. Biringen and A. Saati, “Comparison of several finite-difference meth- 
ods”, J. Aircraft 72, 90 (1990). 


25 


14 G. Erlebacher and M. Y. Hussaini, “Numerical experiments in supersonic 
boundary layer stability”, Phys. Fluids A 2, 99 (1989). 

15 E. M. Saiki, S. Biringen, G. Danabasoglu, and C. L. Streett, “Spatial 
simulation of secondary instability in plane channel flow: comparison of 
k- and h-type disturbances”, J. Fluid Mech. 253 , 485 (1993). 


26 
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Azimuthal 

Modes 

Axial 

Modes 

Linear 

Theory 

p 

0.973203 

0.973227 

0.975964 

pU T 

0.978329 

0.978330 

0.975964 

pue 

0.974753 

0.974760 

0.975964 

pu z 

0.975856 

0.975860 

0.975964 

Ex 

0.973270 

0.973287 

0.975964 


Table 1: Comparison of the growth rates obtained from the DNS and the 
linear theory. Growth rate is the imaginary part of the temporal eigenvalue 
in the linear theory and it is calculated from the time history of the maximum 
amplitudes of the first Fourier modes in the azimuthal and axial directions 
from the DNS. M = 2.0, Re = 2,000, Ri/Ri = 0.5, T 2 /Ti = 1.0, = 

-1.0, a = 19.0. 
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Figure 2: Radial distributions of the velocity, temperature, and pressure 
perturbation fields as obtained from the linear theory corresponding to the 
primary disturbance used in the case studies. Linear-theory prediction for 
the complex frequency is u = 0.5168 + z0.9758. 
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Figure 3: Radial distributions of the velocity, temperature, and pressure 
perturbation fields as obtained from the linear theory corresponding to the 
axial-subharmonic disturbance used in the case studies. Linear-theory pre- 
diction for the complex frequency is w = 0.4517 + i0.7582. 
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Azimuthal Component 



Time 

Figure 4: Initial exponential growth of the maximum amplitudes of the first 
Fourier modes of the flow variables during the direct numerical simulations 
for the harmonic evolution. Thick straight line is the prediction by the linear 
theory. M = 2.0, Re = 2,000, = 0.5, T 2 /T x = 1.0, = 0.0, a = 

19.0, Initial Amplitude = 0.001. 
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Figure 5: Temporal evolution of maximum Fourier amplitudes of harmonics 
in the axial and azimuthal direction. M = 2.0, Re = 2,000, R 2 /R 1 = 0.5, 
and f^/^i = — 1.0. 
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Figure 6: Energy spectra of maximum Fourier amplitudes of harmonics at 
time=25.51. 
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Figure 7: Amplitudes of two-dimensional Fourier modes in the harmonic case 
at different times. 
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Figure 8: Evolution of the mean velocities in the harmonic case for 
t = 12.85, 19.27,25.51 in the direction of the arrows. 
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Figure 9: Contours of Ej at two different cross-sections of (r, z) plane at t = 
0.6313 (a), 12.85 (b), 25.51 (c), 38.00 (d). Contour levels are kept the same 
in (b), (c) and (d). Solid lines are used for positive levels and dotted lines 
for negative levels. Only inner half of the channel is shown for clarity. 
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Figure 10: Isosurfaces of 50% of the maximum of the disturbance total energy, 
Et , at different times during the harmonic case. Only the domain 1.1 < r < 
1.5 is shown for clarity. 
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Figure 11: Temporal evolution of maximum Fourier amplitudes for the sub- 
harmonic case. 
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{N r , No, N z ) = (81 x 130 x 102) 


(iV r , Ng, N z ) = (65 x 66 x 102) 


Figure 12: Comparison of spectral time evolution of the axial-harmonic 
modes up to t = 21 at two different spatial resolutions for the subharmonic 
case. ( N r ,Ne,N z ) = (65 x 66 x 102) case was terminated at t = 21. 
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Figure 13: Radial Fourier amplitude distributions for the mean flow, funda- 
mental mode, and the first subharmonic mode. , t = 15.10, , 

t = 24.95, , t = 29.79. 
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Figure 16: Temporal evolution of the energy transfer rate amplitude (Tin 
Eq. (24) from the fundamental mode (Aq = fc 2 = 1) to the various modes 
(k[,k 2 ) with increasing wavenumbers from top to bottom in each plot. 
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Figure 17: Contours of Et in the subharmonic case at two different cross- 
sections of (r, z) plane at t = 11.00 (a), 21.01 (b), 30.76 (c), 35.64 (d). 
Contour levels are kept the same in all frames. Solid lines are used for 
positive levels and dotted lines for negative levels. Only inner half of the 
channel is shown for clarity. 
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Figure 19: Contour plots of fidug/dr in ( z-0 ) plane on the inner cylinder 
during the subharmonic-development case. 
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Figure 21: Energy spectra of maximum Fourier amplitudes in the forced- 
subharmonic case at time = 20.10. 
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Figure 22: Contours of Ej in the forced-subharmonic case at two different 
cross-sections of (r, z) plane at t = 0.0030 (a), 4.888 (b), 9.844 (c), 14.87 
(d). Contour levels are kept the same in (c) and (d). Solid lines are used 
for positive levels and dotted lines for negative levels. Only inner half of the 
channel is shown for clarity. 
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Figure 23: Contour plots of fidu$/dr in (z-0) plane on the inner cylinder 
during the forced-subharmonic-development case. 
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